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ABSTRACT 

The characteristic Ufetimes of molecular clouds remain uncertain and subject 
to debate, with arguments having recently been advanced in support of short- 
lived clouds, with lifetimes of only a few Myr, and in support of much longer- 
lived clouds, with lifetimes of 10 Myr or more. One argument that has been 
advanced in favour of long cloud lifetimes is the apparent difficulty involved in 
converting sufficient atomic hydrogen to molecular hydrogen within the short 
timescale required by the rapid cloud formation scenario. However, previous 
estimates of the time required for this conversion to occur have not taken into 
account the effects of the supersonic turbulence which is inferred to be present 
in the atomic gas. 

In this paper, we present results from a large set of numerical simulations 
that demonstrate that H2 formation occurs rapidly in turbulent gas. Starting 
with purely atomic hydrogen, large quantities of molecular hydrogen can be pro- 
duced on timescales of 1-2 Myr, given turbulent velocity dispersions and mag- 
netic field strengths consistent with observations. Moreover, as our simulations 
underestimate the effectiveness of H2 self-shielding and dust absorption, we can 
be confident that the molecular fractions that we compute are strong lower lim- 
its on the true values. The formation of large quantities of molecular gas on 
the timescale required by rapid cloud formation models therefore appears to be 
entirely plausible. 

We also investigate the density and temperature distributions of gas in our 
model clouds. We show that the density probabihty distribution function is 
approximately log-normal, with a dispersion that agrees well with the prediction 
of Padoan, Nordlund & Jones (1997). The temperature distribution is similar to 
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that of a polytrope, with an effective polytropic index 7ofr ^ 0.8, although at low 
gas densities, the scatter of the actual gas temperature around this mean value is 
considerable, and the polytropic approximation does not capture the full range 
of behaviour of the gas. 

Subject headings: astrochemistry — molecular processes — ISM: molecules - 
ISM: clouds 



Introduction 



An important goal of star formation research is the development of an understanding 
of the formation of molecular clouds, since these clouds host all observed Galactic star 
formation. The main constituent of these clouds is molecular hydrogen, H2, which forms 
in the interstellar medium (ISM) on the surface of dust grains, a s its formation in the ga s 
phase by radiative association occurs at a negligibly slow rate (iGould fc Salpeten Il963l ). 
The basic physics of II2 formation on dust is believed to be well understood, although 
uncertainties remain about va rious issues such as the relative importan ce of physisorption 
versus chemiso rption (see e.g. iKatz et al\ Il999l : ICazaux fc Tielend I2OOJ) or the role of dust 
grain porosity (jPerets fc Bihanj|2006l ). 



However, many other questions remain unanswered. In particular, there is little consen- 
sus regarding the physical mechanism responsible for cloud formation. One school of thought 
holds that molecular clouds are transient objects, formed and then dispersed on a timescale 
of only a few million years by the action of large-scale turbulent flows in the interstellar 
medium (ISM), which are ther nselves driven by energy i nput from supernovae and from the 
magnetorotational instability (IMac Low fc Klessenll2004j ). Several lines of evidence in favour 
of this picture, such as the absence of po st-T Tauri stars with ages greater than 3 Myr i n 
local star forming regions, are discussed by lHartmann. Ballesteros-Paredes &: Bergiru (120011 ). 



An alternative school of thought holds that clouds are gravitationally bound objects 
in virial equilibrium which evol ve on timescales of the order of 10 Myr or longer, a nd are 
supported by magnetic pressure (ITassis fc Mouschoviasll2004l : lMouschovias et a/.ll2006r). or by 



turbu lence driven by the expansion of internal H ll regions (IMatznei 



2002 



Krumholz et al. 



20061 ). In this picture, i t is ea.sier to reconcile the observed mass of molecular gas in the 



Galaxy (M 
3 Mo yr-i 



El2_ 



IO^Mq; lEvanslll999l ) with the inferred Galactic star formation rate {M^ 



Scald Il986l ) without requiring the average star formation efficiency of molecular 



clouds to be very small. 



One argument often advanced in favour of longer cloud lifetimes is the apparent difficulty 
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involved in producing sufficient H2 in only one or two Myr to explain observed clouds, given 
the relatively slow rate at which forms in the ISM. The H2 formation timescale in the 
ISM is approximately (iHoUenbach. Werner fc Salpeterlll97ll ) 



^form 



10^ yr 



'1^ 



n 



where n is the number density in cm ^, which suggests that in gas w ith a mean numbe r 
density n ~ 100 cm~'^, characteristic of most giant molecular clouds ( IBlitz &: ShdllQSOl ). 



conversion from atomic to molecular form should take at least 10 Myr, longer than the 
entire lifetime of a transient cloud. However, estimates of this kind do not take account of 
dynamical processes such as supersonic turbulence or thermal instability which may exert 
a great deal of influence on the effective H2 formation rate. The impact of these processes 
can be investigated through the use of numerical modelling, but in most models published 
to date either the dynamics of the gas have been highly simplified, with the flow restricted 
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In lGlover fc Mac Low! (120061 . hereafter paper I) we described how we have modified the 
publicly available ZEUS-MP hydrodynamical code (jNormarul2000l ) to allow us to accurately 
simulate the thermal balance of the atomic and molecular gas in the ISM, and to follow 
the growth of the molecular hydrogen fraction with time. An important feature of our 
approach is the use of simple approximations to treat the effects of H2 self-shielding and 
dust shielding, allowing us to approximate the H2 destruction rate without requiring us to 
solve the radiative transfer equation for the photodissociating radiation, which could only 
be done at a prohibitively high computational cost. 

We also presented results from some simple tests of our modified code and showed 
that it produces physically reasonable results. Finally, we used our modified code to model 
molecular cloud formation by the gravitational collapse of quasi-uniform, initially static gas. 
We showed that in this scenario, large-scale conversion of atomic to molecular gas occurs on 
a timescale t > 10 Myr, equivalent to at least 1-2 free-fall times (for simulations performed 
with our canonical value for the initial gas number density, rii = 100 cm~^). 

In the current paper, we present results from a set of simulations performed with a 
supersonically turbulent initial velocity field, and show that the presence of turbulence dra- 
matically reduces the time required to form large quantities of H2. In supersonically turbulent 
gas, the large density compressions created by the turbulence allow H2 to form rapidly, with 
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large molecular fractions being produced after only 1-2 Myr, consistent with the timescale 
required by rapid cloud formation models. We also find evidence that much of the H2 is 
formed in high density gas and then transported to lower densities by the action of the 
turbulence, a phenomenon that may have a significant irn pact on the chemistry of the ISM 
jxie. Allen fc LangeilEoOsI : Iwillacv. Langer fc Allenll2002h . 

The plan of this paper is as follows. In § |2l we briefly describe the numerical method we 
use to simulate the gas and in § E] we discuss the initial conditions used for our simulations. 
In §111 we examine the timescale for H2 formation in turbulent gas, and in §[n]we analyze the 
resulting distribution of the molecular gas. In § |6l we discuss the temperature distribution 
of the gas and its evolution with time, and also examine the behaviour of the effective 
polytropic index, jes- In § [7] we examine how sensitive our results are to variations of our 
initial conditions. We conclude in § [S] with a summary of our results. 



Numerical method 



Pavlovski et al 



We solve the equations of fluid fl ow for the gas in the ISM using a modifled version of 
the ZEUS-MP hydrodynar nical code (NormanI 2000 ). Our modiflcations are based on those 
developed for ZEUS-3D by 
are built on earlier work by 



J2002h and fSmith fc RosenI J2003h . which in turn 



Suttner et all (Il997l ). and involve two major changes to the code. 



First, we have added a limited treatment of non-equilibrium chemistry to the code. We 
treat the chemistry in an operator split fashion. During the source step, we solve the coupled 
set of chemical rate equations for the fluid, together with the portion of the internal energy 
equation dealing with compressional and radiative heating and cooling, under the assumption 
that the other hydrodynamical variables (e.g. density) remain flxed. During the subsequent 
advection step, we advect tracer flelds which track the abundances of the chemical species of 
interest and which advect as densities, while the internal energy is advected as it would be in 
the unmodifled version of ZEUS-MP. For stability, we solve the coupled equations implicitly, 
using Gauss-Seidel iteration. On the rare occasions that the iteration fails to converge, we 
use a more expensive bisection algorithm that is guaranteed to flnd a solution. 

To improve the efficiency of the code, we make use of subcycling. Since the cooling 
time is generally much shorter than the chemical timescale for gas with the range of physical 
conditions of interest in this study, we evolve both the chemical rate equations and the energy 
equation on a cooling timestep, Atcooi, given by: 

Atcooi = 0.3^, (2) 



where A is the net rate at which the gas gains or loses energy due to radiative and chemical 



- 5 - 



heating and cooling: in the convention used in our code, A < corresponds to coohng and 
A > to heating. Often, Atcooi will be m uch shorter than t he hyd rodynamical timestep. 



Athydro, which is calculated as detailed in IStone &: NormanI (Il992a] ). In that case, rather 
than constraining the code to evolve the hydrodynamics on the cooling timestep, we allow 
the chemistry code to subcycle: to take multiple source steps, each of length Atcooi, until 
the total elapsed time is equal to Athydro- The cooling time is recomputed following every 
chemistry and cooling substep, and the code also ensures that the total elapsed time cannot 
exceed Athydro by shortening the last substep (if required) . To further improve the efficiency 
of the code, our implementation of subcycling functions at the level of individual grid zones, 
so that only zones for which Atcooi < Athydro are subcycled. 

In our chemical model, we adopt stand ard solar abundances for hydrogen and helium, 
together with metal abundances taken from lSembach et all (120001 ) . We assume that metals 
with ionization potential less than that of atomic hydrogen (e.g. carbon, silicon) remain 
singly ionized throughout the calculation, and that metals with ionization potentials greater 
than that of hydrogen (e.g. helium, nitrogen) remain neutral. The ionization state of oxygen, 
which has very nearly the same ionization potential as hydrogen, is assumed to track that 
of hydrogen, owing to the influence of rapid charge transfer reactions. We also assume that 
minor ions such as H~, or B.^ can be ignored, and so are left with only four species of 
interest: free electrons, H+, H and H2. The abundances of these four species are constrained 
by two conservation laws - conservation of charge and of the number of hydrogen nuclei 
- and so only two of these species need to be tracked directly. We choose to follow the 
abundances of H2 and H+. The reactions included in our chemical model, and further 
details including justification of our simplifying approximations, are given in Paper I. We 
discuss photodissociation further below. 

The second major modification to the ZEUS-MP code that we have made is to incor- 
porate a detailed treatment of the radiative and chemical heating and cooling of the gas. 
The main heat source in gas with a low level of extinction is photoelectric heating (caused 
by the ejection of photoelectrons from dust grains and polycyclic aromatic hydrocarbons 
illuminated by photons with energies in the range Q < hu < 13.6 eV), while in gas with 
high extinction, c osmic ray heating dor ninates. Our treatrn ent o f these processes is based on 
Bakes fc TieM Jl994h (as modified bv lWolfire era/lboosh and lColdsmith fc Langeil Jl978h 
respectively, with the effects of extinction modelled using the approximation discussed in 
§[2II] below. 

The main coolants at the temperatures and densities of interest in our simulations are 
the fine structure lines of O and Si^. We compute the cooling from these species exactly, 
using atomic data from a variety of sources (see Paper I). We also include heating and cooling 
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from a variety of other processes, which are generally of lesser importance, and which are 
also summarized in Paper I. 



2.1. H2 photodissociation, self-shielding and dust extinction 



Following iDraine fc Bertoldil (Il996l ). we assume that the H2 photodissociation rate can 
be written as the sum of an optically thin rate, fcph,0) and two separate shielding factors, /shield 
and /dust that account for the effects of H2 self-shielding and dust shielding respectively: 



kph — ^ph, /shield /dust- 

For the optically thin rate, we adopt the expression 

fcph,o = 3.3 X 10-^\s-\ 



(3) 



(4) 



where we hav e assumed that the ultraviolet (UV) field has the same spectral shape as the 
Draind (119781 ) field, and where x is a dimens i onless factor which characterizes the i ntensity 
of the field at 1000 A relative to the iHabingl ( 19681 ) field; note that for the original IPraine 
field, X = 1-7. 

To treat the e f fects of self-shielding, we use the approximation for /shield suggested by 
Draine fc Bertoldil Jl996h : 



/ 



0.965 



shield 



+ 



0.035 



[l + x/hy (l + a;)i/2 



■exp [-8.5 X 10-^(1 + a;)^/2] 



(5) 



where x = N^^/b x 10^"^ cm^^, with A^H2 being the H2 column density between the fluid 
element and the source of the UV radiation, and with 65 = 6/10^ cm s~^, where b is the 
Doppler broadening parameter. 

As in paper I, we use two different approximation to compute A^H2 for every grid zone. 
In most of our simulations, we use a local shielding approximation. In this approximation, 
we assume that the dominant contribution to the self-shielding of any given fluid element 
comes from gas in the immediate vicinity of that element, and consequently approximate 
Nu^ as: 

At 

(6) 



Ax 

= -^'^Ha, 



where riHa is the H2 number density in the zone of interest and Ax is the width of the zone, 
measured parallel to one of the coordinate axes. Note that since our grid zones are cubical, 
the choice of axis is immaterial. This approximation is intended to capture, in a crude 
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fashion, the effects on H2 self-shielding of the significant Doppler shifts along any particular 
line of sight. 

In a few simulations, however, we use an approach that we referred to in paper I as 
the six-ray shielding approximation. In this approach, we compute an effective H2 column 
density for each grid zone by first computing exact values along a small number of lines of 
sight and then averaging these values appropriately. To simplify the implementation, the 
lines of sight are chosen to be par allel to the coordinate axes of t he gr id. An approach of th is 



type has previously been used by iNelson fc Langerl (1l997l . Il999l ) and lYoshida et al.l (l2003l ) 



To treat shielding due to dust, we use a similar approach. We follow iDraine fc Bertoldi 



( 119961 ) and write /dust as 

/dust = e"""d.iooo = e"'"d,ioooAfH,tot^ (7) 

where rd,iooo is the optical depth due to dust at a wavelength of lOOOA, o"d,iooo is the corre- 
sponding effective attenuation cross-section (which has a value ad.iooo — 2 x 10~^^ cm^ for 
the diffuse ISM), and A^H,tot is the total column density of hydrogen nuclei between the zone 
of interest and the source of the UV. To compute A^n.tot, we use the same approximations as 
described above. In other words, when we use the local shielding approximation, we assume 
that the dominant contribution to iVn^tot comes from local gas, and so write it as: 

Ax 

A^H,tot = — ('^H+«H+ +2nH2), (8) 

while in the six-ray approximation, we again compute it exactly along a small number of 
lines of sight. 

The main advantage of our local shielding approximation (aside from its simplicity and 
speed) is the fact that it allows us to be certain that we are underestimating the true amount 
of shielding and hence overestimating the H2 photodissociation rate in simulations that use 
it. We can therefore be confident that the H2 fractions computed in these simulation are 
lower limits on the true values and hence that the timescales for H2 formation that we find 
are upper limits on the actual values. As we will see later, this fact serves to strengthen the 
conclusions that we draw from these simulations. 

A major disadvantage of our local shielding approximation is that they make the pho- 
todissociation rate depend explicitly on the numerical resolution Ax. Consequently, the 
equilibrium H2 abundance, XRa.eq, also becomes resolution dependent, as can easily be seen 
from the following equation for XHa.eq 

-n, (9) 



1 •^H2,eq Rph 
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where n is the number density of hydrogen nuclei, -Rfomi is the formation rate of H2 on 
dust grain surfaces, Rp^ = fcph^H2) and where the H2 abundance is defined such that 
= Irin^/n. It must also be said that the physical justification for treating dust shielding 
using the local approximation is somewhat lacking since the shielding provided by the dust 
should not be significantly affected by Doppler shifts along the line of sight. In paper I, 
we saw that for initially static, uniform density gas, these drawbacks are very serious and 
results obtained using the local approximation are of questionable accuracy. However, in 
the turbulent flows studied in the present paper, the performance of the local shielding 
approximation is very much better, as we will see in § H] & [5] below. 

A major disadvantage of the six-ray approximation is the extremely coarse angular 
resolution of the radiation field that it provides. This poor angular resolution will cause us 
to overestimate the amount of shielding in some regions, and underestimate it in others: the 
precise details will depend on the particular form of the density field, but in general we will 
tend to underestimate the amount of shielding whenever the volume filling factor of dense 
gas is small. On the other hand, the fact that this approach does not take account of velocity 
structure along any of the lines of sight means that it may significantly overestimate /shield 
in a supersonically turbulent flow. It is therefore difficult to determine whether the amount 
of H2 produced in simulations using this method is an overestimate or an underestimate 
of the true amount. For the main problem that we are interested in investigating - the 
determination of the H2 formation rate in dynamically evolving, cold atomic gas - this is 
problematic, as it may lead us to derive an artificially short timescale for II2 formation. In 
view of this, we have adopted the local shielding approximation in most of our simulations, 
and have run only a few simulations using the six-ray approximation for the purposes of 
comparison. As we will see later, the results of these simulations agree surprisingly well with 
those of simulations using the local approximation. 

Finally, we note that to compute the effect of the dust shielding on the photoelectric 
heating rate, we again use a very similar approximation: we use our value of A'^n.tot computed 
above to calculate Ay, the extinction of the gas in the V band, and then compute the 
photoelectric heating rate using a radiation field strength attenuated by a factor e"^-^"^^, as 



suggested by iBergin et al\ ( 120041 ) . 
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Initial conditions 



3.1. Box size and initial number density 



Since our aim in this paper is to model the transition from atomic to molecular gas, 
we have chosen to consider relatively small volumes, which we visualize as being smaller 
sub-regions wi thin larger , grav itati onally collapsing clouds o f gas, such as those found in the 
simulations of iKravtsovi (120031 ) or iLi. Mac Low fc KlessenI (120051 ). We therefore consider a 
cubical volume of size L, and apply periodic boundary conditions to all sides of the cube. 
As we discussed in paper I, in choosing a value for L we aimed to strike a reasonable balance 
between simulating a large, representative region of a cloud (which argues for large L) and 
maximizing our physical resolution for any given numerical resolution (which argues for small 
L). For most of the simulations presented here, we settled on L = 20 pc as an appropriate 
value, but in § l7.1l we examine the effects of varying L. The value of L used in each simulation 
is summarized in Tabled! together with the values used for our other adjustable parameters. 

Within the box, we assumed an initially uniform distribution of gas, characterized by 
an initial number density rii. In most of our simulations, we take rti = 100 cm~^ as this 
value is consistent with t he inferred mean densities of many observed molecular clouds 
(IMac Low fc KlessenI |2004| ). However, in § 17.41 we explore the effects of r educing rzj. Atomic 
gas w ith n = 100 cm~^ lies well within the cold neutral medium regime (jWolfire et a/.lll995 



20031 ) and has a short cooling time (tcooi < 0.05 Myr), and so our results are insensitive to 
the initial temperature Ti. In most of our simulations, we adopt a rather arbitrary initial 
temperature T[ = 1000 K, but as we demonstrate in § 16.11 simulations with Ti = 100 K 
produce essentially identical results for times t ^ 0.05 Myr. 



3.2. Magnetic field strength 



Since there is now considerable observati onal evidence for the presence of dyri amically 
significant magnetic fields in interstellar gas (IBeckl 1200 ll : iHeiles fc Crutcherl |2005| ) , we in- 
cluded a magnetic field in the majority of our simulations. For simplicity, in simulations 
where a field was present, we assumed that it was initially uniform and oriented parallel to 
the z-axis of the simulation. The strength of the field was a free parameter, and the values 
used in our various simulations are summarized in Table [TJ Observational determinations 
of the local magnetic field strength give a typical value of 6 ± 2/iG, and so we ensured that 
our fiducial value for the initial magnetic field strength, -Bj^M = 5.85/iG, was consistent with 
this value. 



- 10 - 



Our choice of this shghtly odd value was motivated by our desire that the mass-to-flux 
ratio of the gas should be some simple multiple of the critical value (INakano fc Nakamura 
1978h 

' ^ (10) 

crit 



at which magnetic pressure balances gravity in an isothermal slab. For a ful ly atomic cloud 
of ga s, the mass-to-fiux ratio can be written in units of this critical value as (jCrutcher et al 
2004( 1 

T\/r T\T 

(11) 



— = 3.8 X 10-^^—. 
$ B 



where is column density of atomic hydrogen in units of cm~^ and B is the strength of 
the magnetic field, in units of ^G. For a simulation with a box size of 20 pc and an initial 
atomic hydrogen number density ni = 100 cm^'^, this gives M/$ = 23.45/5^, where Bi is 
the initial magnetic field strength (in units of /iG), and so if Bi = -Bj^M, then (M/ ^)f\f\ = 4. 



Obser vatio ns of magnetic field str engths in molecular cloud cores, summarized in ICrutcher 
(I1999I ) and ICrutcher et al\ (120041 ). find a smaller mean value, (M/$)jnean — 2, but there 
is significant scatter around this value, and in any case it is far from clear that we would 
expect to find the same value of M/$ in dense cloud cores as we would find in the much 
lower density neutral atomic gas. 

In addition to this fiducial case, we ran simulations with a number of other initial field 
strengths, as detailed in Table [1] and in § 17.21 



3.3. Initial velocity field 



We generate the initial v elocity field required f or ou r simulations of turbulent gas by 
using the method described in lMac Low et al\ (119981 ) and lMac Lowl (119991 ). In this method, 
we produce values for each component of velocity in each grid zone by drawing values from a 
Gaussian random field constructed to have a fiat power spectrum for wavenumbers k = L/X^ 
in the range 1 < < 2 and zero power outside of that range. The normalization of the 
resulting velocity field was then adjusted until the root mean square velocity, frms, matched 
a previously specified initial value, frms,i- 

In most of our simulations, we set frms,i = 10kms~^ . Our choice of this v alue is motivated 
by observations of molecular clouds on scales ~ 20pc (jSolomon et a/.lll987l ) which find line- 
of-sight velocity dispersions a\os — 2-7kms~^. If due primarily to turbulent motions in the 
clouds, these velocity dispersions imply RMS turbulent velocities in the range Vrms — 3.5- 
12 km s~^. Our chosen value lies towards the top end of this range, but in § 17.31 we examine 
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the effects of adopting smaller initial values. 

Finally, it should be noted that in the simulations presented in this paper, we consider 
only decaying turbulence, i.e. turbulence which is not maintained by a regular input of 
kinetic energy, a nd which therefore largely dissipates within a few turbulent crossing times 
( iMac Lowlll999t ). However, we find that significant amounts of H2 form within only 1-2 Myr, 
before the RMS velocity has decayed by more than a factor of two. We therefore anticipate 
that the results for driven turbulence would be quite similar. We hope to investigate this 
further in future work. 



4. The H2 formation timescale 

As in paper I, we first consider the formation of H2 in some detail in a set of fiducial 
simulations differing only in numerical resolution, before going on to consider in later sections 
the effects of changing the various input parameters. Parameters for the full set of simulations 
are listed in Table [TJ The values we adopted for our fiducial runs are an initial density 
rti = 100 cm~^, an initial temperature Ti = 1000 K, an initial RMS turbulent velocity 
^^rms,i = 10 km s~^, an initial magnetic field strength Bi = 5.85/iG and a box size L = 20 pc. 
With these parameters, there is initially ~ 0.1 Mj of gas in the box, but this number rapidly 
increases to ~ 10 Mj as the gas cools to a thermal equilibrium temperature of approximately 
65 K. This initial period of cooling takes place in less than 0.05 Myr, just as in the initially 
static simulations discussed in paper I. This is a much shorter timescale than the turbulent 
crossing time of the box, tcross = -^/"^rms.i — 2Myr, and so the gas reaches thermal equilibrium 
before the turbulence has had time to strongly perturb the density structure. The initial 
mass-to-fiux ratio of the gas in the box is M/ $ = 4, in units of the critical value. 

We performed four runs with these parameters, with numerical resolutions of 64^, 128^, 
256^ and 512^ zones, which we designate as MT64, MT128, MT256 and MT512 respectively. 
In Figure [H we plot the evolution with time in these four runs of the mass-weighted mean 
molecular fraction, (a;H2)M5 defined as 

Mm = , (12) 

where pi{2{i,j, k) is the mass density of H2 in zone k), AV{i,j, k) is the volume of zone 
k), Mh is the total mass of hydrogen present in the simulation, and where we sum over 
all grid zones. 

As Figure [1] makes clear, the H2 abundance evolves rapidly in turbulent gas. Large 
quantities of H2 are produced within a short time, with mass-weighted mean molecular 
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fractions of approximately 40% being produced within only 2 Myr. For comparison, in our 
fiducial static runs in paper I, we found that (xh2)m exceeded 40% only for t > 30.8 Myr 
in the run performed using the local shielding approximation (run MS256 in the notation of 
paper I), and at t > 7.6 Myr in the run performed using the six-ray shielding approximation 
(run MS256-RT). The latter time is probably the better estimate, but even in this case there 
is a factor of 3-4 difference in the H2 formation timescale between the static and turbulent 
runs. Indeed, H2 forms significantly faster in our fiducial turbulent runs than in a static 
run in which there was no ultraviolet background and hence no photodissociation of H2 (run 
MS256-nr; in this run, the time at which (xh2)m > 0.4 is again approximately 7.6 Myr). It 
is also apparent from Figure [1] that although the evolution of (a;H2)M in these simulations 
shows some dependence on the numerical resolution of the simulation, for reasons which we 
will explore later, the effect is not large, and the simulations converge on very similar values 
for (xh2)m after the first couple of megayears of evolution. Figure [T] therefore demonstrates 
one of the main results of this paper: the timescale for H2 formation in turbulent gas is much 
shorter than in quiescent gas with the same mean density and large quantities of molecular 
gas can be produced in turbulent regions within only a few megayears. As we will see as we 
proceed, we obtain qualitatively similar results for a wide range of initial conditions. 

Examination of the evolution of the peak density, Pmax, and the RMS density, Prms, in 
our turbulent simulations (plotted in Figure [2] with thick lines and thin lines respectively) 
gives us a strong hint as to why (xh2)m grows so rapidly in turbulent gas. The turbulence 
very quickly produces strong compressions in the initially uniform density field, significantly 
increasing the RMS density of the gas within a megayear, and producing large peak densities. 
At times t > 1 Myr, the growth in prms and pmax slows as the turbulence starts to become 
fully developed, before accelerating again at t ~ 2 Myr (in the high resolution 256^ and 512'^ 
runs) or t ~ 4 Myr (in the lower resolution 64^ and 128^ runs) due to the onset of runaway 
gravitational collapse. 

It is clear from Figure [2] that the values of prms and pmax that we obtain from our 
simulations have not converged, and therefore that the density field is not fully resolved 
even in our 512^ run. This is a major cause of the differences apparent in the evolution 
of (a;H2)M in the various runs - in the higher resolution runs, we resolve more of the dense 
structure, and so form more H2. We would expect this trend to continue as we increase the 
numerical resolution until we reach a point at which the flow is fully resolved. However, as 
discussed in § 4.3.2 of paper I, to fully resolve shocks at the densities of interest requires 
a physical resolution of less than 0.001 pc, which we cannot realistically achieve with a 
fixed grid code. Fortunately, these resolution issues do not seem to have a large impact on 
the values we obtain for (a;H2)M (other than in run MT64, which is very poorly resolved), 
suggesting that our results are not particularly sensitive to the accuracy with which we can 
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follow the densest parts of the flow. As we shall see in § I5.H this is probably because dense 
gas in our simulations quickly becomes almost fully molecular, at which point it clearly 
ceases to contribute significantly to the growth of (xh2)m- Hence, even fairly large errors in 
the computed densities in these fully molecular regions have little impact on (a;H2)M- 

One obvious question is whether the rapid growth in the value of (xh2)m in our simula- 
tions is driven directly by gravitational collapse, which we expect to occur much earlier in a 
supersonically turbulent medium due to the density enhancements created by the turbulence, 
or whether it is these density enhancements alone that are responsible for the increased rate 
of H2 formation. Figure [2] suggests that at early times, it is turbulence which plays the major 
role, with gravitational collapse becoming important only at later times. 

To investigate this point in more detail, we ran a set of simulations in which we disabled 
the effects of self-gravity in the code. The results of these simulations, which we denote as 
MT64-ng, MT128-ng, MT256-ng and MT512-ng are summarized in Table[2], while in Figure[3] 
we compare the evolution of (xh2)m in the highest resolution simulation, MT512-ng, with its 
evolution in run MT512. From the figure it is clear that gravitational collapse is responsible 
for very little of the H2 formation seen in the simulations. Most of the H2 that forms in our 
turbulent simulations does so in dense regions that are not gravitationally bound. This is an 
important point, as it means that the rapid H2 formation that we see in our simulations is not 
specific to the case of gravitationally bound clouds: we would expect to see a similar effect 
in unbound atomic regions in which significant turbulenc e is pre sen t (such as the thin, dense 



sheets formed in the simulations of lAudit fc Hennebelld (120051 ) or IVazquez-Semadeni et al 



( I2OO6I )) provided that enough gas is present to allow for effective self-s hielding of the resultiri g 



H2. Indeed, a similar effect has been noted in a different context by lPavlovski et al\ ( 12002| ). 
who simulate the effects of high-speed (frms > 15 km s~^), decaying turbulence in dense gas 
(n = 10^ cm~^), and show that following an initial period of dissociation, the H2 reforms 
rapidly in the turbulent gas. 

It is also important to establish how sensitive our results are to the approximations 
that we have made in order to treat the effects of the UV radiation field. In Figure UK we 
compare the time evolution of (xh2)m in run MT256 with its evolution in run MT256-RT, a 
run performed using the same initial conditions as run MT256, but which used the six-ray 
shielding approximation instead of the local shielding approximation. It is clear that the 
difference between the two runs is small. Slightly more H2 forms in run MT256-RT, which is 
just as we would expect given the greater amount of shielding in that run, but the difference 
between the two runs is no more than 30% at t ~ 3 Myr. In this plot we also show how 
(a^H2)M evolves in the absence of an ultraviolet background, using results from run MT256-nr, 
which used the same initial conditions as run MT256, save that x = O-O- The results from 
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this run are very similar to those from run MT256-RT. 

This figure demonstrates that for supersonically turbulent flow, our approximations 
work reasonably well. At early times, the rate of growth of (xh2)m is determined primarily 
by the time taken to form the H2, which is essentially the same in all three simulations, and 
so the results agree very well. At later times, H2 photodissociation begins to have more of an 
impact in run MT256, and the results of this run slowly diverge from those of runs MT256- 
RT and MT256-nr. The continued agreement of the latter runs indicates that there is enough 
shielding in run MT256-RT to make (xh2)m almost insensitive to the value of x {^^ least for 
the range of values examined here). It may be that the results from this run are closer to the 
true behaviour than the results of run MT256. However, as we cannot be certain of whether 
the value of (xh2)m we obtain from run MT256-RT is an overestimate or an underestimate, 
whereas we can be certain that the value in run MT256 is an underestimate, we prefer to 
focus on the latter results, since they allow us to draw conclusions that are more robust. 

A further point to investigate is how accurately our simulations can follow (xh2)v, the 
volume-weighted H2 fraction, defined as 

(a^H2)v = , (13) 

where xn^ii.j, k) is the H2 fraction in the grid zone with coordinates (i, j, k), and where N 
is the total number of grid zones. In a turbulent simulation, where the distribution of mass 
is highly inhomogeneous (see § l5.ip . (xh2)v will be far more sensitive to the behaviour of gas 
in low-density regions than (xh2)m- The evolution of (a;H2)v in runs MT256, MT256-RT and 
MT256-nr is plotted in Figure Hb. It is clear from the figure that in this case there is more 
disagreement between the results obtained with our two different shielding approximations, 
and hence more uncertainty in the true value of (xh2)v- However, even in this case, the 
uncertainty is no more than a factor of two at t ~ 3 Myr. The fact that we find better 
agreement for (a;H2)M than for (xh2)v implies that our approximations do a better job of 
modelling the behaviour of high density, H2-rich regions than of low density, H2-poor regions. 
However, this is just what we expect: in regions in which XH2,eq ~ 1, even large errors in /shield 
or /dust will cause only small errors in the H2 fraction, whereas in regions with a;H2,eq <^ 1, 
small errors in /shield or /dust can lead to large errors in the H2 fraction. 

Finally, it is necessary to address the questi on of whether gr a yitati onally unstable regions 



in our turbulent simulations are well-resolved. iTruelove et al\ (119971 ) showed that in order 
to properly resolve collapse and avoid artificial fragmentation, it is necessary to resolve the 
local Jeans length of the gas by at least four grid zones. In other words, collapse is resolved 
only while 

Ax<^Lj(p,T), (14) 
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where Ax is the width of a single grid zone, and Lj is the Jeans length, given by 
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where Cg is the adiabatic sound speed. Since the densest gas in our simulations is also the 
coolest and so has the smallest sound speed, we can determine when (and if) the Truelove 
criterion is first violated by following the evolution with time of Lj in the grid zone with the 
highest gas density. We have done this for each of our simulations, and list the results in 
Table [2l We see that in most of our simulations, the Truelove criterion is violated within the 
first 1-2 Myr. This would appear to call into question our results at later times. However, as 
we discussed previously in paper I, the fact that we no longer properly resolve gravitational 
collapse in dense gas once the Truelove criterion is violated does not necessarily invalidate all 
of our subsequent results: the key question is how much gas is found in unresolved regions. 
To quantify this, we examined intermediate output dumps of density, internal energy and 
H2 fraction from each of our standard runs, determined which zones were unresolved in 
each case, and computed the fraction of the total gas mass in resolved regions, /res, and the 
fraction of the total H2 mass in the same resolved regions, /res,H2 fo^' every output time for 
each run. The values of /res and /res,H2 that we found at the end of each run are summarized 
in Table [HI 

We found that in most of our simulations, we resolved zones containing more than 90% 
of the mass and of the H2 for the whole duration of the simulation, and that in our highest 
resolution run (MT512) we resolved more than 99.8% of the mass and 99.4% of the H2. The 
only run in which our resolution was significantly poorer than this was our lowest resolution 
64^ run MT64, in which we resolve no more than 50-60% of the gas mass and H2 mass. We 
therefore argue that any errors we make in following the subsequent evolution of regions that 
fail to satisfy the Truelove criterion will not have a major effect on the evolution of (xh2)m 
and so will not significantly alter our main conclusions. 



In addition to determining the timescale for H2 formation, we would also like to know 
how the resulting H2 is distributed in the gas. We can get an immediate visual impression of 
this by looking at how varies within a planar slice taken from the simulation volume. For 
example, in Figures [5^ and[5]D we show slices through run MT512 at time t = 1.9 Myr in the 
x-y and x-z planes respectively. We have made use of the periodic boundary conditions to 



5. 



H2 distribution 



5.1. Morphology 
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orient the datacube so that both figures are centered on the zone with the largest density in 
the simulation. These figures demonstrate that the distribution of molecular gas in the box 
is highly inhomogeneous, with the largest concentrations being found in thin, filamentary 
structures that fill only a small fraction of the total volume. Comparison of the two figures 
suggests that there is a higher degree of structure in the x-z plane than in the x-y, possibly 
due to the fact that the magnetic field is initially aligned parallel to the 2;-axis, but the 
highly ordered structure found in the simulations of magnetized collapse without turbulence 
analyzed in paper I is clearly absent. 

An alternative way of visualizing the H2 distribution is to look at the molecular gas 
in projection, by computing the H2 column density distribution that would be seen by a 
distant observer. In Figures and [Hb, we show the H2 column density distribution that 
would be seen by an observer looking along the y axis or the z axis respectively. These images 
again highlight the inhomogeneous nature of the H2 distribution, and also demonstrate that 
considerable small-scale structure exists in this distribution, despite the fact that our initial 
turbulent velocity field included power only on large scales. 

5.2. Density dependence of the H2 fraction 

We can study the H2 distribution in a more quantitative fashion by examining how the 
H2 fraction varies with density. To do this, we computed and n for each of the zones in 
run MT512 at t = 1.9 Myr and then binned the data by number density, using bins of width 
0.05 dex. We then computed the mean and standard deviation for in each bin. The 
resulting values are plotted in Figure [3 Note that although the mean values we compute 
here are volume weighted, we would not expect the mass weighted values to differ greatly, 
since the narrow width of our density bins means that there is little variation in the gas mass 
from zone to zone within a given bin. We also plot in Figure [7] the mean value of XH2,eq for 
each density bin, computed using the local shielding approximation. 

Three features of Figure [7] deserve particular comment. First, it is clear that there is 
considerable scatter in the value of at a given density. Comparison with the plots of 

versus n presented in paper I for smoothly collapsing gas shows that the scatter is much 
greater in the turbulent case. Despite this, however, there is still an obvious underlying 
trend in the distribution of with n, which runs in the direction that we would expect 
(i.e. high density gas is more highly molecular than low density gas). Second, even though 
(a^H2)M — 0.4 at this point in the simulation, there are already regions where the H2 fraction 
is very much higher, and indeed gas with a number density n > 10^ cm~^ is already almost 
entirely molecular. 
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Finally, comparison of the distribution of with n with the equihbrium distribution 
(computed using the same local shielding approximation as is used in the simulation, and 
given by the dashed line in Figure [7]) indicates that in most of the gas, the H2 fraction has yet 
to reach equilibrium. For densities in the range 300 < n < 10^ cm~^, the actual H2 fraction 
lies below the equilibrium curve, indicating that the H2 fraction in this gas is still growing. 
More interestingly, in gas with n < 300 cm~'^, the mean xb_2 hes above the equihbrium curve, 
implying that the H2 fraction in this gas must be decreasing: in other words, more H2 is being 
destroyed by photodissociation than is being formed on dust grains. To confirm that this 
interpretation is correct, we have computed the values of the formation and photodissociation 
rates of H2 for each grid zone in run MT512 at t = 1.9 Myr. The resulting values were binned 
in a similar fashion to the H2 abundance, and the mean values in each density bin are plotted 
in Figure m As we can see, the formation rate is larger than the photodissociation rate in gas 
with n > 300 cm~^, but is smaller than the photodissociation rate in gas with n < 300 cm~^. 

Since our simulations start with fully atomic gas, there are only two ways to explain the 
fact that exceeds XHa.eq at low densities: either XHj.eq must have been much higher in 
the past, or some (or all) of the molecular content of the low density gas must have actually 
formed in higher density gas and then been transported to lower densities by the action of 
the turbulence. Of these two explanations, the first does not appear to be viable - most 
of the parameters that determine the value of a;H2,eq at a given density in our simulations 
(such as the strength of the UV background) are fixed, and the one parameter that does 
vary, the gas temperature, cannot give us the effect we seek. To demonstrate why, we have 
plotted in Figure [H] the variation of XHa.eq with temperature for gas with a number density 
n = 100 cm~^. The mean temperature of the 100 cm~^ gas in our simulation is approximately 
77 K, and the corresponding value of XHa.eq is indicated in Figure [9] by a cross. From the 
figure it is clear that if the temperature were slightly higher at early times in the simulation 
(i.e. if T ~ 200-300 K), then XHa.eq may also have been slightly higher, but at most by about 
20%. On the hand, if the temperature had been much higher, or much lower, then XHa.eq 
would have been smaller than its present value. Given that we find a mean value of x-^^ at 
n = 100 cm~^ that is actually five times greater than the equilibrium value, it is clear that 
temperature variations at earlier times can be responsible for no more than a small portion 
of this excess, and that most or all of it must therefore be due to the transport of H2 from 
high density regions to low density regions. 

To test this hypothesis, we performed several simulations at 256^ resolution in which we 
artifically suppressed H2 formation in gas with a number density smaller than some threshold 
value nth- In runs MT256-th3e2, MT256-thle3 and MT256-th3e3 we used values for nth of 
300, 1000, and 3000 cm~^ respectively. Our standard set of input parameters were used for 
all three runs. In Figure [TOl we show how (xh2)m evolves in these runs; we also plot the 
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corresponding evolutionary track from run MT256 for comparison. Unsurprisingly, we see 
that as we increase rith, the total amount of H2 that forms in the simulations decreases. 
However, the reduction in (xh2)m is less than a factor of two for rith < 1000 cm~^, suggesting 
that most of the H2 that forms (and persists) in the simulation does so in the high density 
gas. 

To examine how the H2 in these simulations is distributed with density, we computed 
the mean value of x^^ in a set of density bins of width 0.05 dex in all three runs at time 
t = 1.9 Myr, and compared the values with those previously computed for run MT256. The 
results are plotted in Figure [TTl for clarity, we do not indicate the size of the standard 
deviations in these values, although they are comparable to those in Figure [3 

There are two important conclusions that we can draw from this figure. First, it is clear 
that although the imposition of a density threshold for H2 formation substantially reduces 
the amount of H2 present in gas with n < nth, it does not eliminate it: some H2 remains, 
and since it cannot have formed in situ in these simulations, it must have formed in higher 
density gas and been transported to the low density gas by the turbulent velocity field. 
Second, imposing a density threshold also affects the H2 abundance at densities above the 
threshold: although the results of the various runs eventually converge at high density, there 
are clear differences for densities n < Srith, rather than at n < nth as we might have initially 
expected. 

In retrospect, this behaviour should perhaps have been expected. Since the H2 formation 
rate scales as the square of the number density, n^, overdense gas makes a disproportion- 
ately large contribution to the overall H2 formation rate compared to underdense gas or 
gas with densities near the mean value. However, previous numerical work on isothermal 
turbulence has shown that many of the overdense structures formed in supersonically tur 



bulent flows are transient objects wi th short lifetimes (see e.g. iKlessen. Heitsch fc Mac Low 



2OOOI : IVazquez-Semadeni et a/.ll2005l ). Therefore, any H2 formed within these transient over- 
densities will not remain at high densities for long, but will quite quickly find itself carried 
into much lower density surroundings. The results of runs MT256-th3e2, MT256-thle3 and 
MT256-th3e3 suggest that a large fraction of the H2 in the gas is produced in this fashion. 

It would clearly be of great interest to determine the extent to which H2 formed in 
transient overdensities is actually mixed into lower density material as the overdensities 
expand and are broken up by the flow; i.e. the amount of turbulent mixing of material which 
occurs. Unfortunately, to determine this with an Eulerian code such as ZEUS-MP, we would 
need to use tracer particles to follow individual fluid elements and these are not currently 
implemented in our version of the code. We intend to revisit this issue in future work. 
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Finally, we have investigated how sensitive our results on the density distribution of 
H2 are to the choice of shielding approximation. In Figure [12] we plot as a function 
of n in runs MT256 and MT256-RT at time t = 1.9 Myr. Above n = 1000 cm-^ the 
results of these two runs are almost indistinguishable. Since, as we have just seen, the 
majority of H2 in these simulations forms at densities n > 1000 cm~^, this provides a simple 
explanation for the insensitivity of (xh2)m to our choice of shielding approximation. At 
lower densities, the difference between runs MT256 and MT256-RT is greater. However, for 
100 < n < 1000 cm~^, the mean value of derived in run MT256-RT lies within the range 
of values found in run MT256, and the two runs only differ significantly at n < 100 cm~'^. 
In run MT256, the H2 fraction at n < 100 cm~^ very quickly falls to almost zero, as there is 
not enough shielding at these densities to maintain a large H2 fraction. In run MT256-RT, 
on the other hand, shielding by the surrounding gas enables H2 to survive in some regions, 
even at these low densities, producing the extended tail seen in the distribution. 



5.3. Density probability distribution function and cumulative mass function 



We can further characterize the H2 distribution produced in our simulations by deter- 
mining at which densities the bulk of the molecular gas resides. We have already seen that 
the highest molecular fractions are found in the densest gas, but in order to quantify how 
much of the total molecular mass is located at high densities, we need to characterize the 
density distribution of the gas. To do this, we have computed the mass-weighted density 
probability function (PDF) for run MT512 at t = 1.9 Myr (plotted in Figure [T3|) . For ease 
of comparison with other values in the literature, we plot the PDF here as a function of 
s = ln(n/ni), where ni = 100 cm~^ is the mean number density of hydrogen nuclei, which 
in our simulations is identically equal to the initial number density of hydrogen nuclei. We 
also plot in Figure [13] a log-normal function of the form 
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where Sm and a are the mean and dispersion of the log-normal distribution, respectively; 
note that these are related by |sm| = ct^/2. To construct the log-normal distribution in 
Figure [13] we have set Sm = 1.32, which is the mean of the actual dataset. Figure [T3 
demonstrates that at densities above the mean (s > 0), the density PDF of run MT512 
is close to a log-normal in form. On the other hand, at densities below the mean, clear 
deviations from a log-normal are apparent: there is a pronounced excess of power at s = 
—2, correspondi ng to n ^ 10 cm~^, and a deficit of p ower at low e r den sities. Previous 



investigations by [Passot &: Vaquez-Semadenil (Il998l ) and IScalo et al\ (Il998l ) of the influence 
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of the polytropic exponent 7 on the shape of the density PDF produced by a turbulent flow 
suggested that a highly turbulent flow with 7 < 1 should produce a power-law density PDF 
at densities above the mean, and that a log- normal PDF is only recovered if 7 = 1; i.e. if 
the gas is isothermal. However, this prediction is based primarily on the results from one- 
dimensional si mulations and the predi c ted b ehav iour is not seen in the three-dimensional 
simulations of iLi. Klessen fc Mac Lowl (120031 ) or iMac Low et al\ (120051 ). Our simulations 
have an effective polytropic exponent 7efr < 1 (see § l6.2l below). but we also see no convincing 
evidence for power-law behaviour at high density, suggesting that in this case the results of 
one-dimensional simulations are also not a good guide to the three-dimensional behaviour. 

For driven isothermal turbulence, the re is general agreement that the resulting den- 
sity PDF should have a log-nor mal form (jPadoan et al\ 119971 : iPassot fc Vaquez-Semadeni 
19981 : iNordlund fc Padoanlll999h . However, there is less agreement on the dispersion of this 
PDF. Based on the results of three-dimensional turbulence simulations, iPadoan et al\ (119971 ) 
predict that the dispersion, cr, is related to the RMS Mach number of the flow, M.^ by 



= In (1 + Q^bM^) . 



(17) 



On th e other hand, based on their one-dimensional simulations, iPassot fc Vaquez-Semadeni 
(119981 ) predict that the relationship between a and M. is actually cr ^ A^. A t t = 1.9 Myr, 
the RMS Mach number in run MT512 is = 5.55, and so according to the iPadoan et al. 
prediction, we should find th at a = 1.471, and hence that = 1.082. On the other hand, the 
Passot &: Vaquez-Semadenil prediction gives a = 5.55 and hence = 15.40. In Figure \T3\ 



we compare these predictions with the actual density PDF in run MT512 at this time. We see 
that although neither prediction fits the data perfectly, the IPadoan et all p rediction comes 
very much closer to doing so than the iPassot fc Vaquez-Semadenil (119981 ) prediction; the 
former merely underestimates the amount of dense gas slightly, while the latter dramatically 
overestimates it. 



Since iPadoan et all simulate driven turbulence in isothermal gas without self-gravity, 
whereas we simulate decaying turbulence in non-isothermal gas with self-gravity, it is perhaps 
not surprising that we obtain a slig htly different result. We have examined whether the 
difference that we find between the IPadoan et all prediction and the true PDF is due to 
the inclusion of self-gravity in our simulations by comparing the density PDF we obtained 
from run MT512 at t = 1.9 Myr with that in run MT512-ng at the same time. We find 
no significant difference, suggesting that self-gravity has so far had little effect on the PDF. 
The difference may rather be due to the fact that we are examining decaying rather than 
driven turbulence, since the highest density portion of the PDF appears to be established at 
an earlier time than the one which we examine here, as can be seen from comparison of the 
density PDF at t = 1.3 Myr (the long-dashed line in Figure [T3l) with that at t = 1.9 Myr. 
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Altern atively, the difference may be due to our softer equation of state. iMac Low et al. 



(I2OO5I ) find a slightly higher level of disagreem ent between th e dispersions of the density 



PDFs in their simulations and the predictions of iPadoan et ali using highly non-isothermal 
turbulent gas to model the larger-scale ISM. Further investigation of this point might be 
worthwhile. Nevertheless, whatever the reason for the residual differences, it is clear that 



the iPadoan et al\ prescription provid es a good description of the density PDF of the gas, 



while the lPassot fc Vaquez-Semadenil prescription does not. 



As far as the distribution of molecular gas is concerned, the main lesson to be drawn 
from Figure [13] is that only a small fraction of the gas is found at densities n > 10"^ cm~^. 
Therefore, even though this material is fully molecular, it does not represent a large fraction 
of the total mass of H2. Instead, most of the H2 is found in regions with number densities 
in the range 10^ < n < 10^ cm~^, as can be seen more easily by examining the behaviour of 
the cumulative mass distribution of H2, plotted in Figure [T5b below. 

Finally, as part of our effort to understand the sensitivity of our results to our treatment 
of UV shielding, we have compared the cumulative mass distribution of H2 in runs MT256 
and MT256-RT, performed using the local shielding approximation and the six-ray shielding 
approximation respectively. The results are shown in Figure [TH We see that although 
slightly more of the total mass of H2 is found at low densities in run MT256-RT than in 
run MT256, the difference is small, suggesting that the results discussed above are not 
significantly affected by our choice of shielding approximation. 



5.4. Morphological evolution 

So far, we have focused primarily on understanding the distribution of H2 in our simu- 
lations at a single moment in time. The particular time that we have chosen, t = 1.9 Myr, 
corresponds to slightly less than a turbulent crossing time at our initial RMS turbulent ve- 
locity, and so represents a point in time at which the turbulence in the simulations is already 
well-developed, but has not yet decayed to insignificance. Moreover, it is also early enough 
that we can be reasonably confident that even in those cases where the runaway gravitational 
collapse of one or more dense clumps is occurring in the cloud, widespread star formation 
has not yet had sufficient time to occur, and therefore that our neglect of star formation and 
the associated feedback processes does not lead us to derive misleading results. 

However, the choice of this particular moment in time remains somewhat arbitrary, and 
it is also not immediately obvious how the II2 morphology of the gas at this point in time is 
related to the H2 morphology at earlier or later times. To investigate this, we have therefore 
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examined the distribution of H2 in run MT512 at several different output times between 
t = 0.6 Myr and t = 2.5 Myr. 

In Figure [T^ . we show how the dependence of on n varies as the gas evolves. 
For clarity, we have omitted any indication in this plot of the scatter around the mean 
values at the different output times, but in each case this is similar to the scatter already 
seen in Figure At t = 0.6 Myr, most of the gas is still atomic: the molecular fraction 
exceeds 50% only at densities n > 6000 cm~^, and at this point in the simulation very 
little gas is found at these densities. Furthermore, a significant quantity of atomic hydrogen 
(~ 20%) remains in this highly dense gas. By t = 1.3 Myr, much more H2 has formed, 
as is evident both from Figure [T5k . and from our previous plot of the time evolution of 
(3^H2)m (Figured]), which shows that it has increased from {xii^)m — 0.07 at t = 0.6 Myr 
to (xh2)m — 0.25 at t = 1.3 Myr, an increase of almost a factor of four. We see from 
Figure [T5k that gas at densities n > lO^cm"^ is now more than 50% molecular, and that gas 
at densities n > 10'^ cm~^ is almost entirely molecular. On the other hand, gas at densities 
n < n\ remains almost entirely atomic. Between these limits, there seems to be a simple 
relationship between the mean fraction and the density, although this relationship is not 
well described by a simple power-law fit. At later output times, the dependence of on n 
remains qualitatively similar to the relationship found at t = 1.3 Myr, but the curve shifts 
to progressively lower densities, indicating that as t increases, the mean fraction found 
at a given gas density also increases. 

We have also computed the cumulative H2 mass distribution for run MT512 for the same 
set of output times. This is plotted in Figure [T5b . We see that at all of the output times 
that we examine, the majority of the H2 is found in the density range 10^ <n< 10^ cm~'^. 
As we saw in the previous subsection, only a small fraction of the gas in our simulations 
has densities n > 10^ cm~^ and so little of the total H2 mass is located at these densities. 
On the other hand, while there is a significant amount of gas which has n < 100 cm~^, the 
molecular fraction in this gas is small, and so its contribution to the total H2 mass is also 
small. Figure [T5b suggests that lower density gas contributes a larger fraction of the total 
H2 mass at later times, consistent with the increase in the mean molecular fraction at fixed 
n noted above: gas closer to the peak of the density PDF is becoming more molecular and 
so is contributing a larger fraction of the mass. 

Finally, the fact that it takes ~ 1 Myr before the RMS density in our highest resolution 
runs stops increasing rapidly, as can be seen in Figure [21 suggests that it takes at least this 
long for the turbulence to become fully developed. This implies that our results at times 
t < 1 Myr will be significantly affected by initial transients, and so may be unrepresentative 
of the behaviour at later times. The results presented in Figure [15] certainly appear to be 
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consistent with this conclusion, as the distribution of H2 at i = 0.6 Myr is quahtatively 
different from the distribution at our three later output times. This suggests that were we 
to allow the turbulence to develop fully in our simulations before beginning to follow the H2 
formation, we would find somewhat different results at early times. Specifically, we would 
expect that in this case the H2 formation timescale would be even shorter than that found 
in the present work. However, it must be stressed that although our current set of initial 
conditions are rather artificial, deliberately suppressing H2 formation until after large density 
inhomogeneities have developed is equally artificial. Since our main aim in this paper is to 
put an upper limit on the H2 formation timescale, we believe that our current procedure is 
justified. 



-24- 




Time (Myr) 

Fig. 1. — Evolution of (xHa)]^ with time in four runs with increasing numerical resolution: 
run MT64 (dotted hne), run MT128 (dot-dashed line), run MT256 (dashed hne) and run 
MT512 (solid line). 
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Time (Myr) 

Fig. 2. — Evolution with time of Pmax and Prms in four runs with increasing numerical 
resolution: run MT64 (thick and thin dotted lines), run MT128 (thick and thin dot-dashed 
hnes), run MT256 (thick and thin dashed lines) and run MT512 (thick and thin sohd lines). 
Note that the time sampling of the data used to compute prms is ten times coarser than that 
used to compute Pmax- 
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Fig. 3. — Evolution of (a;H2)M with time in two 512^ zone runs performed with and without 
self-gravity: run MT512 (solid line) and run MT512-ng (dashed line). 
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Fig. 4. — (a) Evolution of (xHa)]^ with time in three 256^ zone runs, performed using our 
fiducial initial conditions. In run MT256 (solid hne) H2 photodissociation was modelled 
using the local shielding approximation, while in run MT256-RT (dashed line) we used the 
six-ray shielding approximation. Finally, in run MT256-nr (dot-dashed line), the strength 
of the ultraviolet background was set to zero. Note that the six-ray approximation results 
follow the zero radiation results very closely, (b) As (a), but for the evolution of the volume- 
weighted molecular fraction, (xh2)v- 
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Fig. 5. — (a) Slice in the x-y plane through 512^ zone run MT512 at time t = 1.9 Myr 
showing the spatial variation of the H2 fraction. We have used the fact that our simulations 
were performed with periodic boundary conditions to shift the image so that the zone with 
the highest gas density hes at the center of the figure, (b) As (a), but for a slice in the x-z 
plane. 
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Fig. 6. — (a) The H2 column density in 512^ zone run MT512 at time t — 1.9 Myr, computed 
for lines of sight parallel to the y-axis of the simulation, (b) As (a), but for lines of sight 
parallel to the 2;-axis. 
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Fig. 7. — H2 fraction as a function of the number density of the gas (crosses) for 512^ zone 
run MT512 at time t — 1.9 Myr. To compute these values, we binned the data by number 
density, using bins of width 0.05 dex, and computed the mean value of x^^ for each bin. The 
standard deviation in the value of in each bin is also indicated, as is the mean value of 
XB.2,eq for each bin (dotted line). 
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Number density n (cm ^) 

Fig. 8. — Mean H2 formation rate (solid line) and H2 photodissociation rate (dotted line), 
plotted as a function of density, for 512^ zone run MT512 at time t = 1.9 Myr. For n < 
300 cm~^, the photodissociation rate exceeds the formation rate, indicating that more H2 is 
being destroyed in this gas than is being formed in situ. 
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Fig. 9. — Value of XHj.eq as a function of the gas temperature T, computed for gas with 
n — 100 cm~^ which is shielded by a total hydrogen column density A^H,tot = 6.0 x lO^^cm"^. 
Note that this is the same amount of shielding as there would be in a zone with this value of 
n in one of our 512"^, L = 20pc simulations. The sharp transition in a;H2,cq seen at T ~ 250 K 
is due to the sudden onset of self-shielding: at higher temperatures, there is not enough H2 
in the gas, even in the equilibrium state, to make the Lyman- Werner bands optically thick. 
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Fig. 10.— Evolution of (a:H2)M with time in 256^ zone runs MT256 (solid line), MT256- 
tli3e2 (dashed line), MT256-thle3 (dot-dashed hne) and MT256-th3e3 (dotted hne). In runs 
MT256-th3c2, MT256-thle3 and MT256-th3e3, H2 formation was artificially suppressed in 
gas with a number density n < nth = 300, 1000, and 3000 cm~^, respectively. 
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Fig. 11. — H2 fraction as a function of gas number density at t = 1.9 Myr in 256^ zone runs 
MT256 (solid line), MT256-th3e2 (short-dashed hne), MT256-thle3 (dot-dashed line) and 
MT256-th3c3 (long-dashed line). No H2 can form in these runs at densities n < nth, where 
nth = 0, 300, 1000 and 3000 cm~^, respectively, which means that H2 present at n < nth 
must have initially formed in denser gas. For reference, we also indicate the mean value of 
2^H2,eq as a function of number density in run MT256 (dotted line). 
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Fig. 12. — H2 fraction as a function of the number density of the gas for 256^ zone runs 
MT256 (crosses) and MT256-RT (squares) at time t ^ 1.9 Myr. The standard deviation in 
the value of xu^ in each bin is also indicated. 
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Fig. 13. — Mass-weighted density PDF of the gas in 512^ zone run MT512 at times t = 
1.3 Myr (long-dashed hne) and t = 1.9 Myr (sohd hne). Several log-normal functions are also 
plotted: one with the same mean as the t = 1.9 My r data, 5rn = 1 - 32 (sh ort-dashed line), a 
second with a mean Sm = 1.082, corresponding to the Padoan et a/.l(1997) prediction (dotted 
l ine) a nd a third with a mean = 15.40, corresponding to the iPassot fc Vaquez-Semadeni 
(119981 ) prediction (dot-dashed line). 
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Fig. 14. — Cumulative mass distribution of H2 with n in 256^ zone runs using the six-zone 
approximation (MT256-RT; solid line) and the local approximation (MT256; dotted line) at 
time t = 1.9 Myr. 
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Fig. 15. — (a) H2 fraction as a function of number density n in 512^ zone run MT512 at four 
different output times: t — 0.6Myr (dotted line), 1.3Myr (dashed line), 1.9 Myr (dot-dashed 
line) and 2.5 Myr (solid line), (b) As (a), but for the cumulative mass distribution of H2 
with n. 
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6. Thermodynamics 

6.1. Gas temperature: evolution and distribution 

In order to verify that our resuhs are not sensitive to the initial temperature of the gas 
in our simulations, we ran a 256^ simulation (run MT256-T100) in which we set = lOOK 
but kept all of the other input parameters the same as in run MT256. The evolution with 
time of the minimum and maximum gas temperatures, Tmin and Tmax, in this run and in 
run MT256 are plotted in Figure Uni We see that following a short initial period of cooling, 
the behaviour of the two runs is indistinguishable. We have also verified that there is no 
significant difference in the evolution of (xh2)m between the two runs. 

It is also interesting to examine how the gas temperature varies as a function of number 
density in these simulations. To do this, we took the output data from run MT512, binned it 
by number density n, and then computed the mean temperature and the standard deviation 
in the mean for each bin. The resulting values are shown in Figure [T71 Comparison of these 
results with those obtained from run MS256 in paper I demonstrates that we recover the 
same mean temperature in our fiducial turbulent run as in the analogous, initially static 
runs, at least within the range of number densities occupied by gas in the latter. This is 
unsurprising, as the short cooling times lead us to expect that most of the gas in both sets 
of simulations should be in thermal equilibrium, but it is good to verify that this is actually 
the case. 

However, the results from our turbulent runs do differ from those from the initially 
static runs of paper I in two important respects. The first is that gas is found at a very 
much wider range of densities in the turbulent runs, thanks to the strong compressions and 
rarefactions produced by the turbulence. The second is that at low densities (n < 10^ cm~^), 
there is considerable scatter in the temperatures found in the gas, as can be seen from the 
size of the error bars in Figure [TTl this scatter is absent in the static runs. The cause of 
this scatter is the large number of strong shocks present in the turbulent simulations, or, 
more precisely, the large post-shock temperatures, which can be as large as a few thousand 
Kelvin, as can be seen from Figure [161 However, since the post-shock cooling regions are 
generally under-resolved in our simulations, it is likely that this scatter is overstated. 
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Fig. 16. — Evolution of the maximum temperature of the gas, T^a.^, and the minimum 
temperature, Tmin in 256^ zone runs MT256 (thick and thin sohd hues) and MT256-T100 
(thick and thin dashed hues). In run MT256, the initial gas temperature Ti — 1000 K, while 
in run MT256-T100, T; = 100 K. 
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Fig. 17. — Mean gas temperature T plotted as a function of the number density n in 512'^ 
zone run MT512 at time t = 1.9 Myr. The data were binned in a similar fashion to that used 
in the construction of Figure [7] above. The standard deviation in the mean value in each bin 
is indicated wherever it exceeds the symbol size. 
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6.2. Effective polytropic index 

Another thermodynamical quantity of interest is the effective polytropic index of the 
gas, 7efr, defined here as 

d In » , ^ 

7efr = jr-^, 18 
dlnp 

where p is the gas pressure. Most simulations of turbulence within molecular clouds assume 
an isothermal equation of state, in which case 7pfF = 1 at all dens i ties. H owever, simulations 



of turbulent fragmentation performed by lLi. Klessen fc Mac Low! ( 120031 ) for g as with a poly- 



tropic equation of state (in which case 7eff ^ 1 but is independent of p) and by lJappsen et al. 



(120051 ) for gas with a piecewise polytropic equation of state {p (x p^^ for p < pc and p cx p'''^ 
for p > Pc, with 7i 7^ 72) show that the outcome of the fragmentation process is dependent 
on the value of 'jes- Gas with a soft equation of state (low 7efr) fragments more easily than 
gas with a hard equation of state (high 7efr), producing a larger number of fragments with a 
smaller characteristic mass. In view of this, it is interesting to examine how 7efr varies with 
density in our simulations. 

The variation in 7efr with n is plotted in Figure [18] for run MT512 at t = 1.9 Myr. To 
compute the values of 'jeS plotted here, we first computed the mean gas temperature as a 
function of density, as in the previous section, and also the mean molecular fraction as a 
function of density, as in § 15.11 We then used these values to compute the gas pressure as a 
function of density, and from this computed ■jeS through the use of equation [181 We see that 
although the gas in our simulations is not a simple polytrope, since 7eff varies with density, 
its rate of change is small for densities in the range 10 < n < 10^ cm~^, suggesting that 
at these densities a polytropic approximation with 7eff ^ 0.8 may actually be adequate for 
many appl ications. This yalue o f 7eff is in reasonable agreement with the value of 7efr = 0.725 



derived by jjappsen et all (12005! ) for these densities from a synthesis of various observational 



and theoretica l sour ces. However, it is significantly smaller than the values computed by 



Spaans fc Silkl (I2OOOI ) for this range of densities. 



At densities n < 10 cm~'^, we find larger values of 7efr, indicating a stiffening of the 
equation of state. However, at these densities the large scatter present in the gas temperature 
imphes that the gas is not well described by an equation of state that is purely a function 
of density and so the values of 7efr that we have derived for this density regime are of 
questionable validity. 

We also discount the large fluctuations we see in 7eff at high densities {n > 10^ cm^'^). 
These are caused by small fluctuations in the mean temperature, with bin-to-bin variations of 
less than a degree, and our treatment of cooling at these densities is not sufficiently accurate 
for us to trust in the reality of these features. 
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6.3. Influence of the cosmic ray ionization rate 



As pre viously discussed in paper I, m e asurements of the cosmic ray io nization rate in 



diffuse gas ([ McCal l al 



2003 



Lisztl 



2003 



Le Petit. Roueff &: HerbstI l2004f) and dense gas 



( ICaseUi et a/.. ,1998. : iBergin et al.. l999i Ivan der Tak fc van Dishoeckl |2000| ) give results that 
differ by an order of magnitude or more. In most of our simulations, both in paper I and here, 
we adopted an ionization rate ( = 10^^^ s~^, consistent with the value measured in dense 
gas. However, in order to establish the sensitivity of our results to the assumed ionization 
rate, we also performed one 256^ zone simulation with ( = 10"^^ s"^, a value that lies near 
the high end of current determinations. Aside from (, the remaining input parameters of 
this simulation, which we designate MT256-CR, were the same as in run MT256, and so 
comparison of these two runs serves to highlight the effects of varying (. 



We compare the evolution of Tmin and Tmax in runs MT256 and MT256-CR in Figure [ 
We see two main effects: both Tmin and Tmax are increased by about 10-20 K at times 0.02 < 
t < 0.2 Myr, and Tmax (but not Tmin) is increased by a factor of a few at times t > 2 Myr. 
Interestingly, there is little difference between the runs in the interval 0.2 < t < 2 Myr during 
which much of the H2 forms. We would therefore expect to find little difference in the H2 
formation rates in the two simulations, and this is borne out by our results. After 2 Myr 
of evolution, run MT256-CR has formed only 2% more H2 than run MT256, demonstrating 
that the influence of variations in ( is small. 



7. Sensitivity to variations of the input parameters 

As with the runs discussed in paper I, it is interesting to ask what happens as we vary 
some of the more important input parameters. The effects of varying the initial temperature 
of the gas have already been considered in the previous section; the effects of varying the box 
size L, the initial magnetic field strength Bi, the initial turbulent velocity dispersion frms.i 
and the initial density rii are examined in § I7.1H7.4I below. 



7.1. Box size 

To investigate the extent to which the quantity of H2 formed in our simulations and 
the rate at which it forms depend upon our choice of box size L, we have performed further 
256^ zone simulations with L = 10, 30 and 40 pc, which we have designated as runs MT256- 
LIO, MT256-L30 and MT256-L40 respectively. The evolution of (a;H2)M in these simulations 
is plotted in Figure [20b . along with the results from run MT256 for comparison. We see 
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Fig. 18. — Value of 7cfr as a function of n in 512^ zone run MT512 at time t — 1.9 Myr. The 
data were binned as indicated in the text. 
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Fig. 19. — Evolution of the maximum temperature of the gas, Tmax, and the minimum 

temperature, Tmin in 256^ zone runs MT256 (thick and thin solid lines) and MT256-CR (thick 
and thin dashed lines), which were performed with cosmic ray ionization rates C, = 10~^^s~^ 
and C, = 10~^^ s~^ respectively. 
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that as we decrease L, we increase the amount of H2 to form within the first 1-2 Myr of 
the simulation. This is most hkely due to the fact that the network of dense sheets and 
filaments that characterizes the density field in our turbulent runs takes less time to develop 
in simulations run with a smaller L, since the turbulent crossing time of the box is directly 
proportional to L. Evidence for this can be seen in our plot of pmax versus time in Figure [20b. 

In runs MT256, MT256-L30 and MT256-L40 the differences resulting from the change 
in L are small, and have mostly vanished by t = 3 Myr, with the value of (xh2)m in ah three 
runs having converged at about 0.4. The behaviour of run MT256-L10, on the other hand, 
differs markedly: (xh2)m reaches a peak value at t ~ 1.6 Myr of (a;H2)M ~ 0.3 and then 
begins to decline, decreasing by approximately 20% by t = 3 Myr. The reason for this fall 
off in (xh2)m is hinted at in Figure I^Ub : in the L = 10 pc simulation, none of the gas in 
the dense shocked regions undergoes gravitational collapse, and so the peak density begins 
to decrease as the turbulence decays. The total amount of dense gas also decreases as the 
turbulence decays, and so increasing amounts of H2 formed in high density gas finds itself in 
moderate or low density regions where H2 photodissociation occurs faster than H2 formation. 
Consequently, the total amount of H2 present in the simulation declines. However, it should 
be noted that since our treatment of H2 photodissociation leads us to overestimate the H2 
photodissociation rate, it also causes this decline to occur much faster in our simulations 
than would be the case in a real cloud. 



7.2. Initial magnetic field strength 

We ran three 256'^ simulations in which we varied only B-^ in order to determine how 
sensitive our results are to the strength of the magnetic field. In the first of the simulations, 
run HT256, we set B\ = 0.0 and so examined purely hydrodynamical turbulence. In the 
second, run MT256-Bx2, we set Bi = 11.7/iG, twice our standard value, while in the third, 
run MT256-BxlO, we set B^ = 58.5/iG, ten times as large as our standard value. The mass- 
to-flux ratios in these three runs are M/<l> = 00, 2, 0.4 respectively; note that the gas in run 
MT256-BxlO is magnetically subcritical and so will not undergo gravitational collapse. 

In Figure [211 we show how (xh2)m evolves in these three runs and how this compares 
to its evolution in run MT256. Clearly, the growth of (xh2)m is more rapid in the B^ = 0.0 
run than in any of the magnetized runs. This is not unexpected: in our B^ 7^ 0.0 runs, the 
magnetic pressure of the gas helps it to resist compression by the turbulence, and so we 
would expect it to take longer to build up the same amount of high density structure in the 
Bi 7^ 0.0 runs as in the Bi = 0.0 run. Indeed, if we compare how the density PDFs of the four 
runs evolve, we flnd that at early times, there is signiflcantly more dense gas in run HT256 
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Time (Myr) 

Fig. 20. — (a) Evolution of {xi{^)u with time in 256^ zone runs MT256-L10 (solid line), 
MT256 (dashed line), MT256-L30 (dot-dashed line) and MT256-L40 (dotted line), which 
had box sizes L = 10, 20, 30 and 40 pc respectively, (b) As (a), but for the time evolution 

of Pmax- 
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than in any of the magnetised runs, as illustrated in Figure [22b. while by t ~ 2Myr the PDFs 
in the four runs have almost converged, at least for densities near the peak (Figure [22b). 
Similarly, if we compare the evolution of prms in the four runs (Figure [25]) . we see that at 
t < IMyr, prms is between 50% and 100% larger in run MT256 than in any of the magnetized 
runs. 

Another interesting feature of Figure [21] is the fact that in the B[ ^ 0.0 runs, the 
relationship between (xh2)m and B[ is not what we might expect: as we increase Bi, and 
so increase the magnetic pressure of the gas, we find that (xho)m also increase s. This 



behaviour is due to an effect previously noted by [Heitsch. Mac Low &: Klessed (j200ll ) in their 
simulations of isothermal MHD turbulence. They note that as the strength of the magnetic 
field increases, the density enhancements found in their simulations become greater. They 
attribute this to the fact that when the magnetic field is strong, it remains highly ordered, 
causing the magnetic pressure to remain highly anisotropic. Gas can therefore flow along 
the field lines, forming dense shocked layers that are oriented perpendicularly to the mean 
field. On the other hand, when the field is weak it is easily tangled by the turbulent velocity 
field, and so the magnetic pressure becomes far more isotropic, making it harder to form 
highly dense regions. Figures [22] and [22] which show that more dense gas is present at early 
times in runs MT256-Bx2 and MT256-BxlO than in run MT256, suggest that a similar effect 
is at work in our simulations. Further confirmation of this comes from an examination of 
the magnetic energy associated with the x, y and z components of the magnetic field in 
runs MT256 and MT256-BxlO, as plotted in Figure [21 In run MT256-BxlO, the magnetic 
energy associated with the z-component of the field remains orders of magnitude larger 
than the energy associated with either the x or y-components throughout the simulation, 
demonstrating that the field remains highly anisotropic. In run MT256, on the other hand, 
the difference between the three components is much smaller, indicating that the magnetic 
field is far more isotropic. 

Finally, the significant H2 formation that occurs in magnetically subcritical run MT256- 
BxlO, which is not gravitationally unstable, is further confirmation that gravitational insta- 
bility does not cause the rapid H2 formation seen in our simulations, and that turbulent 
compressions dominate even in strongly magnetized gas. 



7.3. Initial velocity dispersion 

To examine the extent to which our results depend upon the strength of the turbulence, 
as quantified by the initial rms velocity, Wrms,i, we have run several simulations with differing 
values of frms.i- In Figure [23, we show how (a;H2)M evolves in runs MT256-vl, MT256-v2.5 
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Fig. 21. — Evolution of (xh2)m with time in 256^ zone runs with different magnetic field 
strengths. We plot results for runs HT256 (solid line), MT256 (dashed line). MT256-Bx2 
(dot-dashed line) and MT256-BxlO (dotted line), which have initial field strengths — 0.0, 
5.85, 11.7 & 58.5 /xG respectively. 
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Fig. 22. — (a) Mass-weighted density PDF at t = 0.63 Myr in 256'^ zone runs HT256 (solid 
line), MT256 (dashed line), MT256-Bx2 (dot-dashed hne) and MT256-BxlO (dotted line), 
performed with different initial magnetic field strengths, (b) As (a), but for t = 1.9 Myr. 
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and MT256-v5, which have frms,i = 1-0, 2.5 and 5.0 km respectively. We see that as we 
decrease the strength of the turbulence, we increase the H2 formation timescale, with the 
time taken to reach (xh2)m ~ 0.4 approximately doubling from t ~ 2 Myr to t ~ 4 Myr as 
we decrease Vrms.i from lOkms"^ to 2.5kms~^. However, the impact of varying frms,i on the 
amount of H2 to form during the simulation is surprisingly small: runs MT256, MT256-v5 
and MT256-v2.5 have all formed roughly the same amount of H2 by the time the simulations 
end, and only in run MT256-vl is the final value of (xh2)m significantly smaller. 

The reason we find faster H2 formation in runs with a larger frms,i becomes clear if we 
examine how the RMS density of the gas varies in these runs. In Figure [261 we plot the time 
evolution of prms in all four runs. We see that as we increase frms,!, we find a faster increase 
in Prms at early times. This behaviour is easily understood on the basis of a simple timescale 
argument: runs with a larger frms,i have a shorter turbulent crossing time, tcross = L/vrms, 
and so in these runs less time is required to produce the highly overdense regions where 
most of the H2 forms. At t ~ 3 Myr in run MT256 and at t ~ 4.5 Myr in runs MT256- 
v5 and MT256-v2.5, we see a sudden increase in prms, consistent with the onset of runaway 
gravitational collapse at one or more locations in the simulation volume, but prior to this the 
evolution of the RMS density in these runs is dominated by turbulence. Since the turbulent 
compressions produce a very similar RMS density in all three cases, it is not surprising that 
we also find similar values for (xh2)m- Note also that the spread in H2 formation timescales 
seen in Figure [251 is comparable to the spread in RMS density growth rates in Figure [26| and 
that the values of (a;H2)M in the three runs converge later than the values of prms because 
H2 formation is not instantaneous, causing the growth in (xh2)m to lag behind the growth 

in Prms- 

Finally, although run MT256-vl appears, on the basis of Figure [25l to behave differently 
from the other runs. Figure [2S] suggests that this is not actually the case: it is simply taking 
significantly longer to produce large overdensities and the associated high H2 fractions in 
this run, which is to be expected given the smaller f rms and hence greater turbulent crossing 
time. 



7.4. Initial density 

In Figure [571 we show how the evolution with time of (xh2)m changes as we reduce 
the initial density of the gas in our simulations. Plotted in the figure are results from two 
runs performed using the local shielding approximation, MT256-nlO and MT256-n30, which 
had initial densities of rii = 10 and 30 cm~^ respectively, and one run performed using the 
six-ray shielding approximation, MT256-RT-nlO, which had an initial density rii = 10 cm~^. 
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To ensure that the mass-to-flux ratio M/ $ remained approximately the same in these runs 
as in our runs with rii = 100 cm~^, we reduced B[ to 1.755 fiG in run MT256-n30 and to 
0.585 fiG in runs MT256-nlO and MT256-RT-nlO. All of the other input parameters had 
the same values as in run MT256. 

It is clear from Figure [271that reducing rii has a dramatic effect on the amount of H2 that 
forms in the simulations, particularly in the runs using the local shielding approximation. 
The decrease in n[ by about a factor of three between runs MT256 and MT256-n30 causes 
a roughly proportionate decrease in (xh2)m at late times. Similarly, the reduction in ui by 
an order of magnitude between runs MT256-RT and MT256-RT-nlO leads to an order of 
magnitude decrease in (xh2)m- The difference between runs MT256 and MT256-nlO is even 
more striking, with an order of magnitude decrease in leading in this case to a reduction 
in (xh2)m by more than a factor of 40. 

The reason for this sensitivity is clear if we look at how the density distribution varies 
as we vary rii. In Figure [28l we plot the mass-weighted density PDF for the suite of runs at 
different densities at t = 1.9 Myr. In each case, s = ln(n/ni), where the mean density n; is 
of course different in each run. As we decrease n[, the peak of the PDF moves closer to the 
mean and the amount of overdense gas decreases. Since the mean density has also decreased, 
the net effect is a sharp drop in the amount of high density gas: in run MT256, at this output 
time, approximately 22% of the total gas mass is in regions with n > 10^ cm~^, while in 
runs MT256-n30 and MT256-nlO, this figure decreases to 2% and 0.05%, respectively. The 
density PDFs in the six-ray runs strongly resemble those in the local approximation runs 
with the same initial density, and so a similar effect occurs in this case. 

Since we already know that most of the H2 that forms in our simulations does so in 
dense gas, it is not surprising that substantially reducing the amount of dense gas available 
has a dramatic effect on (xh2)m- It is clear from a comparison of runs MT256-nlO and 
MT256-RT-nlO in Figure [271 however, that the magnitude of this effect does depend on 
our choice of self-shielding approximation: we see a significantly greater reduction when 
using the local shielding approximation than when using the six-ray approximation. This is 
unsurprising given our previous discussion of the dependence of x^^ on n in our turbulent 
models (see § 15.21) . We have already seen that far more H2 survives at low densities when we 
use the six-ray approximation. This low density H2 represents only a small fraction of the 
total molecular mass present in our Ui = 100 cm~^ simulations, but is far more significant 
in simulations with a lower value of n^. Nevertheless, even in the six-ray run there is a 
large reduction in (xh2)m when rij is reduced, and so although the quantitative details of 
the relationship between rii and (xh2)m depend on the treatment of self-shielding, the basic 
qualitative result does not: conversion of H to H2 is strongly suppressed at < 10 cm~^. 
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Time (Myr) 

Fig. 23. — Evolution of prms with time in 256^ zone runs HT256 (solid line), MT256 (dashed 
hne), MT256-Bx2 (dot-dashed hne) and MT256-BxlO (dotted line), performed with different 
initial magnetic field strengths. 
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Time (Myr) 

Fig. 24. — Evolution with time of the x, y and z components of the magnetic energy in 
256^ runs MT256 (sohd hnes) and MT256-BxlO (dashed hues). The initial magnetic field 
strength in run MT256-BxlO is ten times stronger than in run MT256. In both cases, the 
uppermost line corresponds to the z component of the field. 
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Fig. 25. — Evolution of (xh2)m with time in 256^ zone runs with different initial RMS 
turbulent velocities. We plot results for runs MT256 (solid line), MT256-v5 (dashed line), 
MT256-V2.5 (dot-dashed hue) and MT256-vl (dotted line), which had v,^s,i = 10, 5, 2.5 & 
1 km s~^ respectively. 
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Fig. 26. — Evohition of p^.^s with time in 256^ zone runs MT256 (solid line), MT256-v5 
(dashed line), MT256-v2.5 (dot-dashed line) and MT256-vl (dotted line), performed with 
different initial RMS turbulent velocities. 




Time (Myr) 



Fig. 27. — Evolution of (a;H2)M with time in 256^ zone runs MT256-RT (upper solid 
hnc), MT256 (dashed line), MT256-n30 (dot-dashed line), MT256-RT-nlO (dotted line) and 
MT256-nlO (lower solid line), performed with different initial gas densities. 
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Fig. 28.— (a) Mass-weighted density PDF at t = 1.9 Myr in 256^ zone runs MT256 (solid 
line), MT256-n30 (dotted line) and MT256-nlO (dashed hne), performed with different initial 
gas densities. 
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8. Summary 

The simulations presented in this paper demonstrate that the timescale for H2 formation 
in turbulent gas is much shorter than the corresponding timescale in quiescent gas. For gas 
with a mean number density rti = 100 cm~^, and with a magnetic field strength and RMS 
turbulent velocity consistent with observations, we find that ~ 40% of the initial atomic 
hydrogen can be converted to H2 on a timescale of 1-2 Myr. Moreover, since these results 
are derived from simulations that use a local approximation to treat H2 shielding, they 
represent a lower limit on the amount of H2 formed. (Comparison between a run using a a 
six-ray shielding approximation and one with H2 photodissociation entirely absent suggests 
that this lower limit must lie within about 30% of the true value). The H2 formation 
timescale that we obtain is consistent with that required by models in which molecular 



clouds have short lifetimes, of the order of a turbul ent crossing time (see e.g. lElmegreen 



2OOOI . iHartmann. Ballesteros-Paredes &: BerginI I2OOII ). and so criticism of these models on 



the grounds that they do not provide enough time to form the large quantities of H2 required 
is shown to be groundless. 

By comparing the rate of H2 formation in simulations performed both with and without 
self-gravity, we have shown that the rapid growth of H2 is due to the large density enhance- 
ments created by the turbulent compressions. Self-gravity becomes important only at later 
times, as part of the gas goes into runaway gravitational collapse. We therefore predict that 
H2 should form rapidly in any dense turbulent atomic cloud, regardless of whether or not 
the cloud is gravitationally bound and regardless of whether it is magnetic subcritical or su- 
percritical. Our finding th at forms quickly in a supersonically turbulent flow is consistent 



with the previous work of iPavlovski et al\ (|2002| ). who find a similar effect in simulations of 



high-speed decaying turbulence in n = 10^ cm ^ gas. 

We have investigated the distribution of the H2 formed in our simulations, and have 
shown that most is located at densities between 1 and 100 times our mean number density 
72; = 100 cm~^. Gas denser than lO^cm"^ is fully molecular, but contributes little to the total 
H2 mass as only a small fraction of the gas is found at these densities. On the other hand, 
gas that is less dense than 100 cm~^ contributes little because it is still primarily atomic. 

A surprising result of our simulations is that in low density regions {n < 300 cm~^), 
there is more H2 than we would expect if the gas were in photodissociation equilibrium. 
We have shown that this is explained by the fact that a large fraction of the H2 found at 
low densities actually formed at high densities, in gas with n > 1000 cm~^, and that it 
has subsequently been transported to low densities by the action of the turbulence. This is 
consistent with previous findings that a large proportion of the dense structures produced 
in a supersonically turbulent flow are transient objects which are not gravitationally bound 
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( iKlessen. Heitsch &: Mac Lowll2000l : IVazquez-Semadeni et a/.ll2005l ). Our results show that 
these transient density enhan cements have a profound im pact on the H2 chemistry of the low 
density gas. Previous work by lGarrod et al\ ( l2005l . l2006l ) . using a highly simplified dynamical 
model, suggests that such transient density enhancements will also have a large impact on 
the chemistry of many other tracer species, such as CS. We might also speculate that strong, 
transient density enhancements associated with turbulence i n the diffuse ISM are the caus e 
of the puzzling small area molecular structures discovered by iHeithausenI (120021 . 12004 12006| ) . 



To fully characterize the impact of the turbulence on the chemistry of the gas, it is 



important to determine the amount of mixing which occurs (jXie. Allen &: Langerl Il995 
Willacy. Langer fc Allen! |2002| ). Our current simulations have not allowed us to do this 



since they are strictly Eulerian in nature, but we plan to address this point in future work. 

We have also examined the sensitivity of our results to the choice of several of our main 
input parameters: the initial temperature of the gas, Tj, the box size L, the initial magnetic 
field strength the initial turbulent velocity dispersion frms.i, and the initial density rii. 

We find very little sensitivity to the initial gas temperature. The gas in our simulations 
cools rapidly and reaches thermal equilibrium within only 0.05Myr and at later times displays 
no memory of its initial temperature. Indeed, if we look solely at the mean gas temperature 
at times t > 0.05 Myr, then we find that it behaves rather like a polytrope, with T oc n'^'=ff~^ 
and 7eff — 0.8. However, at low gas densities, the scatter of the actual gas temperature 
around this mean value is considerable, and the polytropic approximation does not capture 
the full range of behaviour of the gas. 

Simulations performed with different values for the box size L do show some differences 
in behaviour, with H2 forming faster when L is small. This appears to be due to the fact 
that altering L alters the turbulent crossing time, which in turn alters the time required for 
large density enhancements to build up within the box. Nevertheless, in every case our main 
result - that H2 forms rapidly, on a timescale of a few megayears - remains the same. 

Varying the strength of the initial magnetic field also leads to some differences in be- 
haviour, although again in every case we find the same basic result. H2 forms more rapidly 
and in greater quantities when = than when Bi 7^ 0, apparently because large overden- 
sities develop faster in the former case than in the latter. Surprisingly, H2 also forms more 
rapidly when the magnetic field is strong and the gas is subcritical than when it is weak 
(but still dynamically significant) and the gas is supercritical. This appears to be due to 
the fact that when the field is strong, it remains highly ordered, causing gas to fiow prefer- 
entially along the field lines and leading to the formation of dense shocked layers that are 
oriented perpendicularly to the mean field. An effect of this kind was previously noted by 
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Heitsch. Mac Low fc Klesseru ( 120011 ) in their simulations of isothermal MHD turbulence. 



Decreasing the strength of the turbulence by decreasing frms,i lengthens the H2 formation 
timescale, as it again takes longer to build up the same amount of dense structure. However, 
H2 formation remains rapid compared to quiescent models (cf. paper I) and mass-weighted 
mean molecular fractions of ~ 40% can be produced in less than 4 Myr as long as v^^^s,! ^ 
2.5 km s~^. Since observations of molecular clouds on scales ~ 20 pc find one-dimension 



velocity dispersions of 2-7 km s ^ (ISolomon et a/.l 119871 ). corresponding to values of v 



3.5-12 km s ^, this implies that H2 formation will be rapid in real clouds. 

The input parameter to which our results are most sensitive is the initial number density 
TT-i. Decreasing rii substantially decreases the amount of dense gas produced in the simula- 
tions, which has the effect of dramatically reducing the production of H2. Based purely on 
our results, we would expect significant H2 formation to occur only in clouds with mean 
densities greater than ~ 10 cm~^, although we caution that this is no better than an order 
of magnitude estimate at this point. 

In summary, we suggest that the key to understanding the formation of molecular 
clouds is understanding the process (or processes) that produce large, dense clouds of tur- 
bulent atomic gas, since our results demonstrate that such clouds will very quickly become 
molecular. Our results also show that the many transient overdense regions that are created 
by the turbulence play a central role in bringing about this rapid rate of H2 formation. 
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Table 1. Input parameters used for each simulation. 



Run 


L(pc) 


Hi (cm ^) 


Ti(K) 




Vrms,i (km S 1) 


Notes 


MT64 


20 


100 


1000 


5.85 


10.0 




MT128 


20 


100 


1000 


5.85 


10.0 




MT256 


20 


100 


1000 


5.85 


10.0 




MT512 


20 


100 


1000 


5.85 


10.0 




MT64-ng 


20 


100 


1000 


5.85 


10.0 


1 


MT128-ng 


20 


100 


1000 


5.85 


10.0 


1 


MT256-ng 


20 


100 


1000 


5.85 


10.0 


1 


MT512-ng 


20 


100 


1000 


5.85 


10.0 


1 


MT256-nr 


20 


100 


1000 


5.85 


10.0 


2 


MT256-RT 


20 


100 


1000 


5.85 


10.0 


3 


MT256-th3e2 


20 


100 


1000 


5.85 


10.0 


4 


MT256-thle3 


20 


100 


1000 


5.85 


10.0 


4 


MT256-th3e3 


20 


100 


1000 


5.85 


10.0 


4 


MT256-T100 


20 


100 


100 


5.85 


10.0 




MT256-L10 


10 


100 


1000 


5.85 


10.0 




MT256-L30 


30 


100 


1000 


5.85 


10.0 




iVi ± ZOU-J-J^U 




1 nn 


1 nnn 

lUUU 


o.oo 


inn 




HT256 


20 


100 


1000 


0.0 


10.0 




MT256-Bx2 


20 


100 


1000 


11.7 


10.0 




MT256-BxlO 


20 


100 


1000 


58.5 


10.0 




MT256-V1 


20 


100 


1000 


5.85 


1.0 




MT256-V2.5 


20 


100 


1000 


5.85 


2.5 




MT256-V5 


20 


100 


1000 


5.85 


5.0 




MT256-nlO 


20 


10 


1000 


0.585 


10.0 




MT256-n30 


20 


30 


1000 


1.755 


10.0 




MT256-RT-nlO 


20 


10 


1000 


0.585 


10.0 


3 



Note. — 1: runs with self-gravity disabled; 2: run with x = 0.0; 3: runs using the 
six-ray shielding approximation; 4: runs with H2 formation switched off in gas with 



n < nth 
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Table 2. tj-es, ^f, and associated values of (xh2)m for all runs in Table [H 



Run tres (Myr) (xH2)M(tres) U (Myr) (xH2)M(tf) 



MT64 


1.17 


0.12 


7.92 


0.80 


MT128 


0.98 


0.15 


7.92 


0.40 


MT256 


1.08 


0.18 


3.01 


0.42 


MT512 


1.49 


0.30 


2.54 


0.42 


MT64-ng 


0.35 


0.04 


7.93 


0.36 


MT128-ng 


1.77 


0.27 


7.92 


0.30 


MT256-ng 


1.78 


0.30 


7.92 


0.23 


MT512-ng 


1.81 


0.35 


2.85 


0.41 


MT256-nr 


0.95 


0.18 


3.48 


0.61 


MT256-RT 


0.98 


0.19 


2.85 


0.55 


MT256-tli3e2 


1.08 


0.15 


4.66 


0.35 


MT256-tlile3 


1.08 


0.11 


3.23 


0.28 


MT256-tli3e3 


1.08 


0.05 


2.06 


0.12 


MT256-T100 


1.08 


0.18 


2.06 


0.35 


MT256-CR 


1.08 


0.19 


2.96 


0.43 


MT256-L10 






3.96 


0.18 


MT256-L30 


0.92 


0.11 


3.26 


0.46 


MT256-L40 


0.66 


0.05 


3.01 


0.41 


HT256 


0.60 


0.25 


1.97 


0.57 


MT256-Bx2 


1.46 


0.33 


7.92 


0.21 


MT256-BxlO 






3.19 


0.43 


MT256-V1 


5.13 


0.21 


7.48 


0.41 


MT256-V2.5 


2.66 


0.22 


5.13 


0.48 


MT256-V5 


1.62 


0.20 


5.61 


0.47 


MT256-nlO 






7.92 


5.9 X 10 


MT256-n30 






3.01 


0.1 


MT256-RT-nlO 






2.76 


0.07 
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Note. — ty-cs is the time at which the Truelove criterion is first 
violated during the course of the run; when no value is given, this 
indicates that the criterion was never violated, tf is the time at which 
the simulation was stopped. 
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Table 3. /res and /res,H2 for all runs in Table [H 



Run 


/res I If j 


p ( -1- \ 
Jres.Haltfj 


MTfi4 




n ^^8 

W . W W(J 


MT1 2S 


n Q74 


n 93^ 


IV J. -L z. i_/ yj 




n 99n 






n 994 


MTfi4-np- 


n QQQ 


1 nnn 

X .www 


MTl 98-np- 

IV J. J. ■'■■'■ft 


1 nnn 

X .www 


1 nnn 

X .www 


MT25fi-n(? 

IVJ. J. VJ ■'■-'■>-) 


1 nnn 

X .www 


1 nnn 

X .www 




1 nnn 

X .www 


1 nnn 

X .www 


MT25fi-nr 


n 883 


n 877 


MT25fi-RT 




n 992 


MT2^fi-th3p2 


n QQ8 


n 993 


MT25fi-th1p3 

IVX X VJ uXXXCO 


n QQ7 


n 988 

W . U\J\J 


MT25fi-th3p3 

IVX X i—iLJ\J uXXfJCO 


n QQQ 


n 992 


MT2^fi-Ti nn 

IVX X ^xJ\J X xww 


n QQ7 


n 992 


MT2^fi-r!R 

IVX X ^tJ\J V^XL 


n QQ7 


n 993 


MT256-L10 


i.nnn 


I.nnn 


MT256-L30 


n.987 


n.97i 


MT256-L40 


n.985 


n.963 


HT256 


n.992 


n.986 


MT256-Bx2 


I.nnn 


I.nnn 


MT256-BxlO 


i.nnn 


I.nnn 


MT256-V1 


n.936 


n.847 


MT256-V2.5 


n.968 


n.93i 


MT256-V5 


n.979 


n.953 


MT256-nlO 


I.nnn 


I.nnn 


MT256-n30 


I.nnn 


I.nnn 


MT256-nRT-nlO 


I.nnn 


I.nnn 
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Note. — /res is the fraction of the total 
gas mass that is situated in zones that 
satisfy the Truelove criterion, computed 
at the end of each simulation; /res,H2 is 
the fraction of the total H2 mass situated 
in these zones 



